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Abstract 

A space-time adaptive method is presented for the reactive Euler equations de¬ 
scribing chemically reacting gas flow where a two species model is used for the chem¬ 
istry. The governing equations are discretized with a finite volume method and dy¬ 
namic space adaptivity is introduced using multiresolution analysis. A time splitting 
method of Strang is applied to be able to consider stiff problems while keeping the 
method explicit. For time adaptivity an improved Runge-Kutta-Fehlberg scheme is 
used. Applications deal with detonation problems in one and two space dimensions. 
A comparison of the adaptive scheme with reference computations on a regular grid 
allow to assess the accuracy and the computational efficiency, in terms of CPU time 
and memory requirements. 


1 Introduction 

Real world industrial or environmental problems, e.g., management of industrial risks, typ¬ 
ically involve physical and chemical phenomena having a multitude of dynamically active 
spatial and temporal scales. Their direct numerical modelling thus leads to prohibitive 
computational cost. Introducing adaptivity can be understood in the sense that the com¬ 
putational effort is concentrated at locations and time instants where it is necessary to 
ensure a given numerical accuracy, while efforts may be signihcantly reduced elsewhere. 
Adaptive methods are in many cases more competitive than schemes on regular hne grids, 
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in particular for solutions of nonlinear PDEs exhibiting a non-uniformly distributed regu¬ 
larity of the solution. Reliable error estimators of the solution are essential ingredients of 
fully adaptive schemes. They are based for example on Richardson ideas of extrapolation, 
adjoint problems or gradient based approaches. For evolutionary problems, a major task 
is the time evolution of the grid and its reliable prediction for the next time step. How¬ 
ever, to become efficient, adaptive methods require a signihcant effort on implementing 
data structures, which are typically based on graded trees, hash-tables or multi-domains. 
Moreover, the computational cost per cell is signihcantly increased with respect to uniform 
discretizations. Hence, an adaptive method is only efficient when the data compression 
is large enough to compensate the additional computational cost per cell. Fortunately, 
for problems exhibiting local discontinuities or steep gradients, adaptive computations are 
faster than hne grid computations. 

Adaptive discretization methods for solving nonlinear PDEs have a long tradition and 
can be tracked back to the late seventies [ 12 ]. Adaptive hnite element methods have a long 
history, in particular for elliptic problems. For chemically reactive flow with detailled chem¬ 
ical reaction in three dimensions mm proposed stabilized hnite elements with adaptive 
mesh rehnement. The equations are treated fully coupled with a Newton solver and the 
solution of large linear non-symmetric, indehnite systems becomes necessary for which a 
parallel multigrid solver is used. Moving grid techniques have been applied successfully to 
combustion problems [39] . A posteriori error estimators have also been studied for a long 
time to improve the grid, since the early work of Babuska and Rheinboldt [2]. However, 
adjoint problems have to be solved which are linear although the original PDF can be 
nonlinear [3]. Fully adaptive hnite element discretizations of reaction-dihusion problems 
encountered in electrocardiology have been proposed in [121 E3]- For time adaptivity a 
stepsize control with linearly implicit time integrators is used. In space a multilevel hnite 
element method is combined with a posteriori local error estimators. 

The main challenge is to estimate and control the error of adaptive schemes with respect 
to the exact solution, or at least with respect to the same numerical scheme on an under¬ 
lying uniform grid. Self adaptive methods are preferred as they automatically adjust to 
the solution. The block-structured adaptive mesh rehnement technique (AMR or SAMR) 
for hyperbolic partial diherential equations has been pioneered by Berger and Oliger [7|. 
While the hrst approach utilized rotated rehnement grids that required complicated con¬ 
servative interpolation operations, AMR denotes today especially the simplihed variant of 
Berger and Collela |B] that allows only rehnement patches aligned to the coarse grid mesh. 
The striking efficiency of this algorithm, in particular for 3D instationary supersonic gas 
dynamics problems, has been demonstrated by Berger et al. in |1]. 

Recently, multiresolution (MR) techniques have become popular for hyperbolic conser¬ 
vation laws, going back to the seminal work of Harten [36| in the context of hnite volume 
schemes and cell-average MR analysis. Starting point is a hnite volume scheme for hyper¬ 
bolic conservation laws on a regular grid. Subsequently a discrete multiresolution analysis 
is used to avoid expensive hux computations in smooth regions, hrst without reducing 
memory requirements, e.g. for ID hyperbolic conservation laws (Harten [30]), ID con¬ 
servation laws with viscosity (Bihari [5), 2D hyperbolic conservations laws (Bihari and 
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Harten M), 2D compressible Euler equations (Cliiavassa and Donat [H]), 2D hyperbolic 
conservation laws with curvilinear patches (Dahmen et al [19]) and unstructured meshes 
(Abgrall and Harten [T], Cohen et al US]). A fully adaptive version, still in the context 
of ID and 2D hyperbolic conservation laws, has been developed to reduce also memory 
requirements (Gottschlich-Miiller and Muller [35|, Kaibara and Gomes [iQ], Cohen et al 
|17j). This algorithm has been extended to the 3D case and to parabolic PDEs (Roussel 
et al jlH], Roussel and Schneider M), and more recently to self-adaptive global and lo¬ 
cal time-steppings (Muller and Stiriba jlil, Domingues et al [ 23112311231 ) • Therewith the 
solution is represented and computed on a dynamically evolving automatically adapted 
space-time grid. Different strategies have been proposed to evaluate the flux without re¬ 
quiring a full knowledge of hne grid cell-average values. Applications to shock waves in 
compressible flows addressing the issue of shock resolution have been presented in [5], and 
extensions to the Navier-Stokes equations in the weakly compressible regime can be found 
in [33]. Adaptive MR methods with operator splitting have been proposed for multiscale 
reaction fronts with stiff source terms in [2HI EZj. The numerical analysis of the above 
higher order operator splitting techniques has been performed in [201129] and it was shown 
that the splitting time step can be even larger than the times scales involved in the PDEs. 
Adaptive MR computations using the above method with complex chemistry and includ¬ 
ing detailed transport can be found in [21], further applications involving various stiffness 
levels have been presented in [30] [3T] . 

The MR approach has also been used in other contexts. For instance, the Sparse 
Point Representation (SPR) method was the hrst fully adaptive MR scheme, introduced 
by Holmstrom [38] in the context of hnite differences and point-value MR analysis, leading 
to both GPU time and memory reduction. In the SPR method, the wavelet coefficients 
are used as regularity indicators to create locally refined grids, on which the numerical 
solution is represented and the finite difference discretization of the operators is performed. 
Applications of the SPR method have been published in [ 221133 ] • Discontinuous Galerkin 
methods, they have been applied to hyperbolic conservation laws in [13] using Haar wavelet 
indicators to decide where to refine or to coarsen the meshes. These publications reveal that 
the multiresolution concept has been applied by several groups with success to different 
stiff problems. For comprehensive literature about the subject, we refer to the books of 
Gohen [15] and Muller [33] . 

The objective of the paper is the extension of the adaptive multiresolution method [38l 
123] to the numerical simulation of detonation waves. As model we use here a one-step 
chemical reaction involving two chemical species only. Since the chemical source term is 
stiff, a time splitting has to be made between the convective and source terms and each 
term needs to be computed with a different time step and a different time integration 
method. An error-controled time step based on a Runge-Kutta-Fehlberg approach is used 
for the source term integration. The application concerns Ghapman-Jouguet detonations 
in one dimension with different stiffness values and instabilities of detonations waves due to 
an interaction with a pocket of partially burnt gases in two dimensions. These detonation 
problems motivate the use of the Euler equations, although extensions of the method to 
compute viscous flows using the Navier-Stokes equations is straightforward, as proposed 
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The outline of the paper is the following: hrst, we present the set of reactive Euler equa¬ 
tions for a simplihed detonation model. Then we describe the hnite volume discretization. 
A Strang splitting technique is utilized to account for temporal scales in the source term 
that do not influence the hydrodynamics. We also briefly summarize the multiresolution 
strategy. Finally, we show the numerical results for the test problems in one and two space 
dimensions and we set the conclusions together with some perspectives for future work. 

2 Governing equations 

For modelling the combustion process, we use the reactive Euler equations, as described in 
[I8l [3^ . The simplest description of a chemically reacting gas flow assumes that the gas 
mixture is made only of two chemical species, the burnt gas, denoted with subscript b and 
the unburnt gas, denoted with subscript u. The unburnt gas is converted to burnt gas via 
a single irreversible reaction. We represent the mixture state by a single scalar variable Z 
corresponding the mass fraction of the unburnt gas. We also assume gases in the mixture 
to be ideal polytropic gases with equal specihc heat ratio 7 and specihc gas constant r. The 
system of equations in two dimensions, which has been non-dimensionalized in a suitable 
way, may be written as 


dQ dF dG 
dt dx dy 
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Here p denotes the mixture density, V = {vx,Vy)^ the mixture velocity, e the mixture total 
energy per unit of mass, p the pressure, T the temperature and k the chemical reaction 
rate. The two equations of state completing the model are 
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where Qo denotes the amount of heat per unit of mass released in the chemical reaction. 

The reaction rate k{T) of the irreversible chemical reaction is expressed in Arrhenius 
form as 



(5) 


where the pre-exponential coefficient A and the activation temperature T 4 are empirical 
constants. When the reaction source term is stiff, however, the reaction rate may be 
simplihed by adopting the so-called ignition temperature kinetic model, i. e. 



T>Ti 
T <Ti 


( 6 ) 


where Tj denotes the ignition temperature and r the characteristic time of the chemical 
reaction, which determines the stiffness of the problem. This formulation has been chosen 
in the applications presented in the numerical results section. However, the numerical 
method is not limited to this simplihed model and its extension to the Arrhenius law or 
even more complex chemical reactions is possible. 

3 Numerical method 

In this section, we hrst describe the classical Strang splitting and the space discretization 
of the convective terms. Subsequently, the time integration is discussed and, hrst, for 
the convective terms, a time step depending on the CFL is chosen, then, for the source 
term, we split the convective time step into error-controlled time steps using an explicit 
Dormand-Prince method. A diherent choice has been proposed by Duarte et al in [28] . 
based on implicit and explicit Runge-Kutta methods and a posteriori error estimators. We 
also briehy recall the multiresolution method, previously published in ^ 81 123] . 

3.1 Strang splitting 

We denote by C{Q) the operator of the convective terms. Equation Q becomes 


-^ = C{Q) + S{Q). 


(7) 


Discretizing explicitly with hrst order in time, we get 

QU+l + ^(Q")] 


( 8 ) 


where n denotes the time instant and At the convective time step. 
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The splitting relies on the separation of the convective and source terms operators. To 
get a hrst-order discretization, it writes 

Q* = Q^ + AtS{Q^) (9) 

Q-+1 = Q* + AtC'(Q*) (10) 

g- ( 11 ) 

where denotes the hrst-order accurate source term operator with time step At and 
the hrst-order accurate convection operator with time step At. 

To get second-order accuracy, which corresponds to the classical Strang, the following 
procedure can be applied 

/An+l _ qo 2 /^o2 qo 2 /-i 

V — ‘^At/2 ^At ‘^At/2 V 

where S°^^i 2 denotes the second-order accurate source term operator with time step At/2 
and the second-order accurate convection operator with time step At. 

We note that the splitting time step in the above method is hxed. Techniques to 
introduce adaptive splitting time steps based on a posteriori error estimators have been 
introduced in [201 EZ] allowing automatic error control. 


3.2 Space discretization of the convective terms 

Discretization in space is made using a hnite volume method, to ensure conservative hux 
computations. Convective terms are discretized using the AUSM-I- scheme [12]. In this 
procedure, pressure terms are computed separately. The Euler hux F writes 
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where M denotes the Mach number, c the speed of sound and h the enthalpy per unit of 
mass. Denoting by $ the purely convective term and by If the pressure term, we discretize 
the Euler hux F in space following 
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The interface speed of sound is Cj_,_i = y'CjCj+i and the interface convective term is 


$ 


i+i 


if Mj+i > 0 
otherwise 


( 15 ) 


The terms and Pj_,_i follow 
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where 
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(19) 


Third-order accuracy of the spatial discretization far from discontinuities is obtained 
using a MUSCL interpolation EH, together with a Koren slope limiter BH. Denoting by 
q one of the conservative quantities (he. density, momentum, energy, partial mass of the 
unburnt gas), the corrected value via a MUSCL interpolation is 

(li = (li + ^0 (u) {qi - Qi-i) + ^0 {qi+i - qi) (20) 


and 


q'i+i = qi+i - -(j) (u+i) {qt+i - qf 



{qi +2 - qi+i) 


( 21 ) 


where 


ri = 


qi+i - qi 


qi - qi-i 


( 22 ) 


and 

/ l + 2r 
0 , min f 2r, —-—, 2 

Nevertheless, since we use a second-order accurate Strang splitting in time, the global 
accuracy of the scheme remains second-order. 



0 (r) = max 
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3.3 Time integration of the convective terms 

Time integration of the convective term is made using a classical third-order TVD Runge- 
Kutta scheme, i.e. 



= Q” + AtC'(Q”) 

(24) 

Q** 

= ^ [3 Q" + Q* + Af C-(Q*)] 

(25) 

jn+l 

= J [Q" + 2 g** + 2 At (F (g**)] 

o 

(26) 


The corresponding Runge-Kutta tableau is given in Table 
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Table 1: Butcher tableau corresponding to the compact TVD third-order Runge-Kutta 
method. 


The convective time step At is chosen to satisfy the Courant-Friedrichs-Levy (CFL) 
condition. In the computations, CFL = 0.5 is used. 

3.4 Time integration of the stiff source terms 

Due to the stiffness of the chemical source terms in case of detonation computations, a 
high order time integration is chosen to ensure the numerical stability. As the chemical 
time scale is much smaller than the convective one we apply a much smaller time step 
Atd than for the convection terms. However, this is only required in a small area of the 
computational domain, i.e. in the reaction zone. In the rest of the computational domain, 
the source term is almost equal to zero. 

For this reason, we introduce an adaptive local time step Atd — ^ being the 

number of substeps required by the source term computation. In order to adapt this local 
time step with time, a fourth-hfth order embedded Runge-Kutta method is used. The 
computational error between the fourth-order and the hfth-order method enables to decide 
wether the local time step needs to be increased or decreased, as proposed by Fehlberg 
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In this article, an embedded explicit Dormand-Prince method [26] was chosen. The coeffi¬ 
cients are computed in order to minimize the error of the hfth-order solution. This is the 
main difference with the Fehlberg method, which was constructed so that the fourth-order 
solution has a small error. For this reason, the Dormand-Prince method is more suitable 
when the higher-order solution is used to continue the time integration. The coefficients 
are given in Table 
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Table 2: Butcher tableau corresponding to the Dormand-Prince method. The hrst row of 
b coefficients gives the fourth-order accurate solution, and the second row yields order five. 
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The local relative error of the quantity q is denoted by e and the accepted tolerance by 
The optimal time step a Atd can be determined by multiplying the scalar a times the 
current time step At^. The scalar a is given by 


^ ^ / e Atd Y 

\ 2 e At-fYiax / 

In practice, however, the new time step At'^ 


^ /e At^Y 
~ \ e At) 

is determined in function of a, following 



r 2 Atd a a > 2 

= l ^ if^<l 

I Atd otherwise 


(32) 


3.5 Multiresolution method 

The principle in the multiresolution (MR) setting is to represent a set of function cell 
averages as values on a coarser grid plus a series of differences at different levels of nested 
grids. The differences contain the information of the function when going from a coarse to 
a hner grid. 

A tree data structure is an efficient way to store the reduced MR dataas it allows 
to reduce the memory with respect to a hnite volume (FV) scheme on the hnest level. 
This representation could also increase the speed-up during the time evolution because it 
reduces the time-searching of the elements. 

In the following we consider a hierarchy of regular grids in 2D, Qe, 0 < i < L. The root 
cell is Do, 0,0 = D and corresponds to a rectangle with side lengths and hy. The different 
node cells at a level £ > 0 forming are denoted by where {i,j) E A^. The ensemble 
of indices of the existing node cells on the level £ is A^. Note that are rectangles 

with side lengths hx,£ = 2~^hx and hy/ = 2~^hy. In the tree terminology, the rehnement 
of a parent node cell at level £ produces four children nodes D£+i, 2 i, 2 j+i, 

D£+i, 2 j+i, 2 j and D£+i, 2 j+i, 2 j+i at level £ -|- 1, as illustrated in Figure]^ 


^l+l.(2i.2j+li / ^t22M,2J+ll 



Figure 1: Dyadic grid rehnement in 2D. 
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The cell-average value of the quantity u on the cell is given by . j u(x, y) dx dy. 

and correspondingly the ensemble of the existing cell-average values at level i hj Ui = 
{'d£,i,j)(i,j)eAr The projection (or restriction) operator 


■ Ui+l I—>■ U£. 


estimates the cell-averages of a level i from the ones of the level I + 1. The parent cell- 
average is the weighted average of the children cell-averages 

^£,1,] = ^(h£+l,2i,2j + U£j^i^2i,2j+1 + 'h£+l,2i+l,2jr- + 'Uf+l,2i+l,2j-|-l) 

and thus the projection operator is exact and unique. To predict the cell-averages of a 
level ^ +1 from the ones of the level we use a prediction (or interpolation) operator 

Pt^l+l '■ U£ H->• U£+l. 


This operator yields an approximation f/^+i of at the level I + 1. In this paper, we 
use third order interpolation given by a tensor product approach [5]. For n,p E {0,1}, we 
dehne 


^^+l,2i+n,2j+p 


o ( Ij) 

O 


+ 8 ^ 1 )^ U£^ij-i) 

T( 1) ^ [("^£,1+1 j+i ('U^^j—ij-i-i 'U£^j_ij_i)]. 


(33) 


First, this prediction is local, since it is made from the cell average U£^ij and the eight 
nearest uncles Second, it is consistent with the projection, i.e. P£+i^£o= 

Id. The difference between the exact and the predicted values at three children cells yields 
the wavelet (or detail) coefficients. The sum of the four details in the children cells is equal 
to zero [9] 


d£+l,2i,2j+l = 'd£+l,2i,2j+l ~ 'Wr+l,2i,2j+l 

d£+l,2i+l,2j = h£-|-l,2i+l,2jr- — Mr+l,2i+l,2j (34) 

d£+l,2i+l,2j+l = 'd£+l^2i+l,2j+l — U£^i^2i+l,2j+l- 

This implies that the knowledge of the cell-average values on the children is equiv¬ 
alent to the knowledge of the cell-average values on the parents U£ and the wavelet coef- 
hcients De+i = {d£+i^ 2 i, 2 j+i,d£+i, 2 i+i, 2 j,d£+i, 2 i+i, 2 j+i){i,j)&Ae- The so-called multiresolution 
transform on the cell-average values is obtained by repeating this operation recursively on 
L levels [36] . 

Ul <—t {Dl, Dl_i, ... ,Di, Uq). 
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In conclusion, knowing the cell-average values of all the leaves t/n is equivalent to knowing 
the cell-average value of the root Uq and the details of all the other nodes of the tree 
structure. 

In the MR scheme, instead of using the representation on the full uniform grid the 
numerical solution U^r = is formed by cell averages on an adaptive sparse grid 

pn _ r^. Grid adaptivity in the MR scheme is related with an incomplete tree structure, 
where cell rehnement may be interrupted at intermediate scale levels. This means that T” 
is formed by leaf cells 0 < i < L, {i,j) G which are cells without children. 

Here C{Ai) denotes the ensemble of indices for the existing leaf cells of the level i. 

Three basic steps are undertaken to evolve the solution from to 

Refinement: 

-n+l - 

Evolntion: EmrU^j^ 

Coarsening: ^ T{e)UM^ 

The rehnement operator R is a precautionary measure to account for possible translation 
or creation of hner scales in the solution between two subsequent time steps. As the 
regions of smoothness or irregularities of the solution may change with time, the grid T” 
may not be convenient anymore at the next time step Hence, before doing the time 
evolution, the representation of the solution should be extended onto a grid T”'"'', which is 
expected to be a rehnement of T”, and to contain Then, the time evolution operator 

Emr = EMR^At) is applied. The subscript MR in Emr means that only the cell-averages 
on the leaves of the computational grid T"’’*' are evolved in time, and that an adaptive hux 
computation Emr{U^ji) is adopted at interfaces of cells of diherent scale levels. Finally, a 
thresholding operation T(e) (coarsening) is applied in order to unrehne those cells in 
that are unnecessary for an accurate representation of The choice of the threshold 

value is motivated by an error analysis equilibrating the perturbation and discretization 
errors. For details we refer to [miiH]. 

To compress data in an adaptive tree structure, while still being able to navigate 
through it, gradedness is required. For instance, for a given node in the dynamic tree 
structure we assume that: 

• its parent and eight nearest uncles are in the tree (if not, we create them as nodes); 

• for flux computations, if is a leaf, its four nearest cousins ^le,i± 2 ,j and ^le,i,j ±2 in 
each direction are in the tree (if not, we create them as virtual leaves); 

• if a child is created, all its brothers are also created; 

For more details of these procedures, we refer to |1H] . 

In the tree structure, the thresholding operator T(e) is dehned by removing leaves 
where details are smaller than a prescribed tolerance e, while preserving the gradedness 
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property, and the refinement operation R adds one more level as secnrity zone, in order 
to forecast the evolntion of the solntion in the tree representation at the next time step. 
These two operations are performed by the following procednre. 

We denote by A the ensemble of indices of the existing tree nodes in T”’*', by /^(A) the 
restriction of A onto the leaves, and by A^ the restriction of A to a level £, 0 < £ < L. For 
the whole tree, from the leaves to the root: 

• Compnte the details on the nodes {i,j) G A^-i, by the mnltiresolntion trans¬ 

form; 

• Define the deletable cells, if the details on the corresponding nodes and their brothers 
are smaller than the prescribed tolerance. 

For the whole tree, from the leaves to the root: 

• If a node and its children nodes are deletable, and the children nodes are simple 
leaves (withont virtnal children), then delete their children. 

• If the node and its parents are not deletable, and it is not at the maximnm level, 
then create the children for this node. 

To illnstrate the adaptive finx compntation, we consider the leaf f2£+i,2i+i,2j, sharing an 
interface with another leaf at a lower scale level, as illnstrated in Figure For 

the calculation of the outgoing numerical flux on the right interface, we use the cell width 
in the x direction hx/+i as step size. The required right neighboring stencils are obtained 
from the cousins ^£+ 1 ^ 21 + 2 , 2 j and Vt£^i^ 2 i+z, 2 j+ii which are virtual cells. For conservation, 
the ingoing flux on the leaf is set equal to the sum of the outgoing fluxes on the 

neighbour leaves of level ^ + 1. For more details on the implementation of this procedure 
we refer to |48j . 



Figure 2: Adaptive numerical flux computation in 2D. 
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4 Numerical results 


In all numerical simulations, we consider an initially planar Chapman-Jouguet detonation 
moving with constant unitary speed through the unburnt gas to the right of the domain. 
Setting 7 = 1.4, r = 1, Qo = 1 and Ti = 0.22 and assigning the following burnt gas values 


Pb 'y ) ^x,b ^ 1 Pb 1 ) ^b 0 ) 


(36) 


where 


h = 


/ 2(7-1) 

7 + 1 


(36) 


the corresponding von Neumann state past the shock wave is [SniET] 


Pn — 


7 


1 -h 


, = 26 , Pk = 1 + 7(5 , Zk = 1 , 


(37) 


and the unburnt gas values are 


7 


Pu 1 I r 5 '^X,U 0 ) Pu 1 7 ^ 5 1 • 

1+0 


(38) 


The resulting temperature of the unburnt gas = 0.215995 is only slightly lower than 
the chosen ignition temperature Tj. The initial condition for an initial front at x = Xq is 


(' 3 . Q') ^ I - qb) exp [a{x - Xo)] +qB if x < xq 

^ \ qu ii X> Xq 


(39) 


where a is set to K Here q stands for p, v^, p, and Z, while Vy = Q everywhere. This way 
the initial condition is close to the classical Zeldovich-von Neumann-Dring (ZND) solution 
of the reaction Euler equations. 


4.1 One-dimensional detonation 

First we consider a one-dimensional setting. The computational domain is H = [—3,1]. 
In this case, the location of the initial front is set to xq = —2. In order to follow the 
detonation wave and perform longer computations, a change of variables is made for Vx- 
We set v'^ = Vx — S and omit the superscript ' everywhere. Since the initial transition 
happens in the segment [—3, —1], we only focus on the domain [—1,1]. The dimensionless 
hnal time is t = 2.2. 
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4.1.1 Non-stiff case 

In this first part, we consider a detonation with a time coefficient r = 10“^ . As observed in 
Figure the detonation front propagates from left to right with constant maximal values 
of density, pressure and velocity, which correspond to the von Neumann state. The slopes 
for the density, pressure, and velocity, are moderate on the left side of the shock, i.e. in 
the burnt gas zone, but become sharp on the right side, at the interface with the unburnt 
gases. For all the displayed quantities, the curves £t with the reference computation, which 
has been performed on a hne grid with L = 14 scales, i. e. with 16384 grid points. 

Table assembles results for computations performed with different maximum number 
of scales for both the MR and the FV method. The relative error is computed on the 
von Neumann value of the density. The numerical value is averaged in time to damp its 
oscillations and is compared to the theoretical one. We observe for the MR computations 
a CPU time compression rate growing with increasing number of scales, while the relative 
error remains comparable with the one of the same computation on the regular hne grid. 
As expected, the error is reduced when increasing the number of scales L. 


Method 

L 

e 

CPU time 

CPU compression 

relative 

error 

MR 

9 

10 “^ 

1 s 

41,34 % 

1.099- 

10-3 

FV 

9 


3 s 

100,00 % 

1.109- 

10-3 

MR 

10 

5■10-3 

4 s 

34,28 % 

6.519- 

10-2 

FV 

10 


13 s 

100,00 % 

6.601 - 

10-2 

MR 

11 

2.5 ■ 10-3 

11 s 

21,86 % 

3.554 - 

10-2 

FV 

11 


54 s 

100,00 % 

3.617- 

10-2 

MR 

12 

1.25 ■ 10-3 

40 s 

18,75 % 

1.740 - 

10-2 

FV 

12 


3 min 37 s 

100,00 % 

1.783 - 

10-2 


Table 3: Chapman-Jouguet detonation front: CPU time, CPU compression and relative 
error on the von Neumann density for different scales and for the FV and MR methods, 
r = IQ-h 

4.1.2 Influence of the stiffness 

In this part, we perform computations of stiff problems and set r first to 10“^, then to 
10“^. In order to resolve correctly the detonation front, 13 scales are required for the case 
r = 10“^. The results are shown in Figures and We observe a much sharper and 
thinner reaction zone, reduced to a few computational points in the case r = 10“^. We 
observe a good agreement with the reference computations, which were computed on a 
regular hne grid with L = 15 scales. 

Table 1^ gives the results of the computations for r = 10“^ and r = 10“^. As expected, 
a high compression rate is reached when 13 scales are used. But as it is seen in Figure 
the error on the von Neumann density remains quite large, even if the detonation front is 
well tracked. 
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Figure 3: Chapman-Jouguet detonation front: density {top, left), pressure {top, right), 
temperature {middle, left), partial mass of the limiting reactant {middle, right), velocity 
{bottom, left) and adaptive mesh {bottom, right) at t = 2.2, L = 11 scales, e = 2.5 ■ 10“^, 
T = IQ-h 
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Figure 4: Chapman-Jouguet detonation front: density {top, left), pressure {top, right), 
temperature {middle, left), partial mass of the limiting reactant {middle, right), velocity 
{bottom, left) and adaptive mesh {bottom, right) at f = 2.2, L = 12 scales, e = 1.25 ■ 10“^, 
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reference - 
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Figure 5: Chapman-Jouguet detonation front: density {top, left), pressure {top, right), 
temperature {middle, left), partial mass of the limiting reactant {middle, right), velocity 
{bottom, left) and adaptive mesh {bottom, right) at f = 2.2, L = 13 scales, e = 6.25 ■ 10“"^, 








Method 

T 

L 

e 

CPU time 

CPU comp. 

rel. error 

MR 

10-2 

12 

1.25 ■ 10-3 

43 s 

20,08 % 

1.336-10-3 

FV 

o 

1 

to 

12 


3 min 37 s 

100,00 % 

1.342 ■ 10-3 

MR 

10-3 

13 

6.25 ■ 10-^ 

2 min 03 s 

7,56 % 

3.057-10-3 

FV 

10-3 

13 


27 min 13 s 

100,00 % 

3.057-10-3 


Table 4: Chapman-Jouguet detonation front: CPU time, CPU compression and relative 
error on the von Nenmann density for different values of the stiffness coefficient and for 
the FV and MR methods. 

4.2 Two-dimensional detonation 

Now we consider a two-dimensional problem. The interaction between a detonation wave 
and a pocket of unburnt gas is simulated. The computational domain is [0,4] x [—1,1], the 
initial planar detonation front is located at x = 1.7 and the initial spot of unburnt gas is 
centered in (xo,|/o) = (2,0). The radius of the circular pocket is tq = 0.1. The parameters 
of the detonation front are the same as in the non-stiff case, i.e. we choose r = 10“^. The 
dimensionless hnal time is t = 1.3. 

In Figure]^ we observe the destabilization of the planar front by the pocket of unburnt 
gas. We remark circular structures due to the expansion of the pocket in the three directions 
and their interaction with the detonation wave going from left to right. Recirculations are 
observed in the center of the pocket, which are advected by the flow. We also observe 
reflexions of the waves on the free-slip boundaries in y = — 1 and y = 1. 

On the right side of Figure we observe that the adaptive mesh follows well the 
structures of the different waves that interact in the flow. A zoom in the center of the 
mesh (Figure enables to see how well the mesh adapts to all the structures of the flow. 
The computation, performed on L = 10 scales i.e. a maximum of 2^^ = 1,048,576 points, 
requires 2 h 59 min of CPU time on a 16 Core PC. The same computation on the regular 
hne grid lasts 14 h 35 min. Hence the CPU time compression is around 20 %. 

Figure shows a cut of the density at y = 0. We observe an excellent agreement 
between the MR and FV computations with L = 10 scales, and a good agreement of both 
curves with the FV reference computation obtained with L = 11 scales, which validates 
the grid convergence of the computation. The total number of grid points in the FV 
computations is 2^^. 

5 Conclusion and perspectives 

We presented an extension of the adaptive multiresolution method for the reactive Euler 
equations, able to deal with fast chemical reactions. A hnite volume scheme with second 
order shock capturing schemes is used for space discretization. Classical Strang splitting 
is applied to deal with the stiffness of the physical problem in time. Space and time 
adaptivity are then introduced using multiresolution analysis and Runge-Kutta-Fehlberg 
schemes, respectively. The implementation uses dynamic memory allocation with tree data 
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Figure 6: Time evolution of the interaction between a detonation wave an a pocket of 
unburnt gas computed with the adaptive MR method. Numerical Schlieren of the density 
gradient {left) and corresponding meshes {right) at f = 0 {top), t = 0.65 {middle) and 
t = 1.3 {bottom). 
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Figure 7: Adaptive MR computation of the interaction between a detonation wave an a 
pocket of unburnt gas: Zoom on the adaptive mesh at t = 1.3. 
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Figure 8: Interaction between a detonation wave an a pocket of unburnt gas. One¬ 
dimensional cuts of density for y = 0 at t = 1.3. Comparison of the adaptive MR and the 
FV computations for L = 10 with the FV reference computation for L = 11. The total 
number of grid points in the FV computations is 2^^. 






































































































































































































































































































































































































structures. 

Applications to detonation problems in one and two space dimensions validated the al¬ 
gorithm and illustrated the efficiency of the adaption strategy. The adaptive computations 
yield both speed-up of CPU time and memory compression with respect to uniform grid 
computations while the precision is automatically controlled using suitable thresholding 
techniques. 

As perspective we plan the extension of the method to three space dimensions and 
also to include multi-species chemical reactions and viscous effects. The parallelization 
of the adaptive method using tree data structures to handle the adaptive grid remains a 
challenging task for the near future. 
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